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Abstract 


Second order turbulence models of the Mellor and Yamada type have been widely 
used to simulate the PBL. It is hov/ever known that these models have several deficien- 
cies, For example, they all predict a critical Richardson number which is about four 
times smaller than the Large Eddy Simulation (LES) data, they are unable to match 
the surface data, and they predict a boundary layer height lower than expected. 

In the present model, we show* that these difficulties are all overcome by a single new 
physical input: the use of the most complete expression for both the pressure-velocity 
and the pressure-temperature correlations presently available. Each of the new terms 
represents a physical process that was not accounted for by previous models. The 
new model is presen :ed in three different levels according to Mellor and Yamada’s 
terminology, with new, ready-to-use expressions for the turbulent moments. We show 
that the new model reproduces several experimental and LES data better than previous 
models. As far as the PBL is concerned, we show that the model reproduces both the 
Kansas data as analyzed by Businger et ah in the context of Monin-Obukhov similarity 
theory for smaller Richardson numbers, as well as the LES and laboratory data up to 
Richardson numbers of order unity. We also show 7 that the model yields a higher PBL 
height than the previous models. 
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1. Introduction 


Reynolds stress turbulence modeling began in the early 40’s (Chou 1940, 1945) and 
since then it has been developed by both physicists and engineers, e.g., Rotta (1951), 
Lumley and Khajeh-Nouri (1974), Launder et al. (1975), Pope (1975), Zeman and 
Lumley (1979), Speziale (1991), and Shih and Shabbir (1992). The parameterizations 
of the turbulence closures have been formulated theoretically, verified experimentally 
(including comparison with the ever more reliable LES data), and applied to various 
engineering flows. In the geophysical applications, Mellor and Yamada (1974, 1982) 
pioneered the use of turbulence closure models to study the Planetary Boundary Layer 
(PBL). The Mellor- Yamada (MY) model and its numerous variants have been more 
successful in the simulation of the PBL than many of the empirical models and have 
been widely used to describe the atmospheric PBL and the oceanic mixed layer. The 
MY models are, however, not without deficiencies. Comparison of MY model results 
with measured data and LES data show consistent discrepancies, and close examina- 
tion indicates that the weakness of the model comes from three sources: (1) a crude 
parameterization for the pressure- velocity and the pressure-temperature correlations. 
(2) the use of a single ’’master” length scale (a'l the length scales corresponding to 
different processes are assumed to be proportional to a master scale), and (3) a down- 
gradient approximation for the third order turbulent moments. These three aspects 
can be handled as three independent components in the model development and each 
of them deserves a separate discussion. Along with many other efforts, the present 
authors tried to address items (2) and (3) elsewhere (Cheng and Canuto 1994; Canuto 
et al. 1994, 2000). The present paper concentrates on item (1), i.e., how to improve 
the parameterization for the pressure correlations, thus generalizing the MY models 
and improving the comparison with both measured and LES data. 

Let’s look at the deficiency (1) of the MY models and its variants (e.g., Hassid and 
Galperin 1983) more closely. Firstly, these models predict too low a critical Richard- 
son number (around 0.2), beyond which the turbulence ceases to exist, while both 
measurements and LES data (e.g., Webster 1964; Wang et al. 1996) indicate that 
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the critical value is around unity. Second, when applied to the surface of the neutral 
boundary layer, none of these models is capable of differentiating between the vertical 
and lateral components of the turbulent kinetic energy, w 2 and v 2 , in fact, they yield 
identical expressions for the two, while experiments consistently show that the vertical 
component is much smaller than the lateral one (Table 1 of Mellor and Yamada, 1982; 
Nieuwstadt 1984, 1985). 

As we will show below, these deficiencies are associated with the oversimplifica- 
tion of the parameterizations of the pressure-velocity correlation fl tJ and pressure- 
temperature correlation Ilf which will be corrected by adopting a more complete ex- 
pression. Both II,j and flf have been shown to contain a slow (return-to-isotropy) part 
and a rapid part (Launder et al. 1975; Lumley 1978) . The rapid parts of both fl tj and 
nf contain velocity terms related to the mean strain-rate tensor S^j, and the vorticity 
tensor R,j. as well as buoyancy terms related to the heat fluxes. In addition, the rapid 
part of flf also contains a term related to the temperature variance 0 2 . By contrast, 
MY models of n tJ include only the slow part and some of the rapid part (the term 
proportional to eS tJ , where e is the turbulent kinetic energy); for flf, only the slow part 
is included. Since each of these missing terms represents a specific physical process, it 
seems appropriate and necessary to incorporate them in the model formulation, as we 
do in the present paper. 

In Section 1, we introduce the general problem. In Sections 2 and 3, we describe 
the basic equations and the new turbulence closure. In Section 4, we derive the general 
model equations for the 3D case. In Section 5, we give the detailed and ready-to-use 
equations for the PBL. The new model is presented in three different ’’levels” according 
to MY’s terminology. Model constants are determined in Section 6. In Section 7, we 
compare the new model and the h4Y model with measured and LES data, where we 
can see that the new model matches the measured and LES data better than previous 
models. Conclusions are presented in Section 8. 
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2. Basic equations 


To model a PBL, we need both mean and turbulent variables. They are: 
a. First moments 
1) Mean velocity, U,: 

DU, d 1 dP n n TT 

Tij - g, - - 2 t ijk UjU k 


Dt dxj 
2) Mean potential temperature, 0: 


p dx. 


where 
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Here u l is the ith component of the turbulent velocity fluctuation, g t = (0,0, g) is the 
gravitational acceleration, P is the mean pressure, p is the mean density, fij is the 
rotation of the Earth, r,y are the Reynolds stresses and hi is the heat flux. 
b. Second moments 
1) Reynolds stresses, r t y: 


D ^ f dUj dU t x 
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Here, a is the volume expansion coefficient, IT t y is the pressure- velocity correlation 
tensor, u is the molecular viscosity and e is the dissipation rate of the turbulent kinetic 
energy e, and D tJ is the diffusion term. 

Of special interest is the equation for the turbulent kinetic energy e: 

1 , 


C= 2 9 


q 2 = U{U{ 


(2e) 
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2) Heat flux, /?,: 


where 


De 1 _ dUi 

Dt + 2°“ ~ TlJ d Xj + P' k ' ' 
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n ? = D < = gir u ' u i e 


dx x ' 
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where Ilf is the pressure- temperature correlation, and D x is the diffusion of the heat 


flux 6,. 

3) Temperature variance, 6 2 : 


7 — 0 2 + D$ = —2h ,—- — 2tg 
Dt ox, 


(4a) 


where 



d 


D e = —u,0 2 

OXi 


(46) 


where \ is the molecular conductivity and D e is the diffusion of the temperature 


variance. 

In the present study, terms containing the molecular viscosity v and molecular 
conductivity \ have been neglected, except for and tg. In addition, in the second 
moment equations, rotation has also been neglected. The modeling of the third-order 
moments exceeds the scope of the present paper, but the interested readers may refer 
to recent work on the subject (Canuto et ah 1994, 2000). As already stated, in this 
paper we concentrate on the closure parameterization of the correlations U tJ and nf, 
which will be shown to improve the PBL model results. 


3. Turbulence closure 


a. Pressure correlations 

The pressure correlation terms n,j and 11^ in Eqs.(2b) and (3b) contain three dis- 
tinct contributions due to (1) turbulence self-interactions (the return-to-isotropy or 
slow part), (2) mean shear-turbulence interactions (a. rapid part), and (3) buoyancy- 
turbulence interactions (also a rapid part). The most complete models for n.j and nf 
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are given by (Launder et al. 1975; Zeman and Lumley 1979): 


Un = n[j J + nj 2) + u\f, n? = n? (1) + n? (2) + nf 


0(3) 


(5a) 


where 
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where b tJ is the traceless Reynolds stress tensor defined as follows 

2e 

b{j — muj v{j 


The other tensors are defined as follows: 
1J 2\dxj dii) 


R ^i_(du 1 _du i 

tJ 2 dx, 


(56) 
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where S tJ and R tJ are shear and vorticity respectively. In most past second-order 
turbulence models for the PBL, these pressure correlations terms were parameterized 
much less completely. For example, the MY model and its variants (Mellor 1973, Mellor 
and Yamada 1974, 1982 and Hassid and Galperin 1983) only consider: 


rd 1 ) _ 2 t -1 b n 2) — --eS n 3) = 0 

tip — ZT P v 0 ui 11 u 5° ,3 ' > 11,J 


n? ,; 1 - rfhi, n'<" = n' |3 > = o (6) 

In other words, only the slow terms and one single rapid term (the first term in the 


t0(2) tt0( 3 ) _ 


expression of n[ 2) ) are retained; most of the rapid terms are neglected, and no buoyancy 
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effects on the pressure correlation are included. In Section 6, we will show how these 
missing terms lead to some of the model deficiencies, e.g., failure to match the data in 
the neutral surface layer as well as in the stably stratified flows. 

4. Algebraic Reynolds stress and heat flux models (AM) 


a. Prognostic equations 

The mean wind U{ and mean potential temperature © are solved prognostically 
using Eqs.(la-c). 

b. A hierarchy of turbulence models 


We present our new AM model as a hierarchy of different levels (levels 3, 2.5 and 2). 
In the level 3 model, the turbulent kinetic energy e and the temperature variance 9 2 are 
solved from their prognostic equations (see Appendix), while the other second moments 
are solved from algebraic equations. In the level 2.5 model, the prognostic model 
equation for 6 2 is reduced to an algebraic relation. In the level 2 model, the prognostic 
equations for both e and 6 2 are replaced by corresponding algebraic equations. 

The t urbulent kinetic energy dissipation rate e and temperature variance dissipation 
rate to could be solved from their prognostic equations, or they could be parameterized. 
Here we express them in terms of the corresponding dissipation time scales, r and tq: 


_ 2e _ 02 

T ’ C ° To 


(?) 


In Section 6 we will discuss a parameterization of r and tq. 

In the main text of this paper we will concentrate on the level 2.5 and 2 models. 


c. Algebraic equations for the second moments 

Combining (2a) and (2f), we can obtain the equation for b tJ , Eq.(5d): 

jjf X Dij = — -eS,j Sjj Z{j + B{j n,j 
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where 


„ d 
1J ~ dx k 


(u t Uj - u h 


(8b) 


Assuming that the left side of (8a) can be neglected and employing (5b) for the 
pressure-velocity correlation Ilq, one obtains the following algebraic equation for 6 tJ 


where 


bjj Aj ct S j j A2 r jj A3T ^ jj | A ^ x ^3 tj 


l= T -f, A 2 = i(l-a,)A 

A 3 = )(1 - a 2 )A, A< = iftA 


(9a) 


(96) 


These model constants will be given in Section 6. Similarly, in the prognostic equation 
(3a) for the heat flux 6 t , if the left side is neglected and (5c) is employed for the 
pressure-temperature correlation Ilf, one obtains the algebraic equation for 6, at level 
3: 

( 2e \ <90 

Ai 3 h 3 = -r (bij + — + Xot/3,0 2 (10a) 

where 

A t j — A5 <5,j + Xqt Sij + A7 t Rjj ( 106) 

T 3 5 

A 0 = 1 — 71, A 5 = — , Ae = 1 — 7<>3? A7 = 1 — -Q3 (10c) 

V 4 4 

At levels 2.5 and 2, we further simplify the problem by neglecting the left side in 
the prognostic equation for d 2 , Eq.(4a), to obtain the algebraic equation 


0 2 = -r e h t ^~ 
ox. 


( 11 ) 


Substituting (11) into (10a), we obtain the algebraic equation for hi at levels 2.5-2: 

2e , \ <90 


w 


here 


and where 


, (, 2e \ <90 

A tJ hj - r (b i} + 3 6 t] j 


<90 

A'ij = X 5 S tJ + A 6 t S tJ + A7 t R ij + AgT 2 /?,-^-- 
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(12 a) 


(126) 


(12c) 


8 


5. PBL Model 


a. Mean wind and potential temperature equations 


In the PBL, several approximations can be made to the equations for the mean 
wind and temperature. In Eq.(la), the horizontal pressure gradient can be expressed 
in terms of the mean geostrophic wind components U g and V g as follows: 


1 

P 


aP <ir)=M v i' -v,) 


dx dy 

and the rotation term can be approximated as 


(13a) 


2tijk£ljUk — 


(136) 


where x, y and z are the eastward, northward and vertical directions respectively, 
f c = 2f \sin<p is the Coriolis parameter with ft the angular velocity of the Earth and 
<? the latitude. In Eq.(lb), the horizontal temperature gradient can be approximated 
with the thermal wind relation, 


d<d dQ\ 


) = A 

( dV 9 

' 

du g \ 

/ 9° 


dz ) 


(13c) 


The equations for the eastward and northward horizontal mean wind components U 
and V and for the mean potential temperature 0 in the PBL can then be written as: 


dQ 

dt 


dU 

dt 


= MV-V g )~ 


duw 
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dV 
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-MU - u g ) - 


dvw 

~d7 


-Tz w0 + 


ga 




(13d) 

(13c) 

(13/) 


b. Level 3 model 

Since the level 2.5 and level 2 models catch the main features of the second-order 
closure models and are easy to use, they have become the most popular second-order 
closure models in the PBL community. We will concentrate on them in the sections 
below. Yet, the level 3 model has its own strength in that it produces counter gradient 
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heat fluxes, a phenomenon observed in the upper part of the PBL. In the Appendix we 
will present the details of the level 3 model for completeness and for future reference. 

c. Level 2.5 model 


In the level 2.5 model, the turbulent kinetic energy e is solved form its prognostic 
equation: 


dt <9 1- 


dU dV 


— = — — — u 2 w -f v 2 w + w 3 — — uw — — vw 4 - gawd — e (14) 

Ot Oz 2 oz oz 


The differential equation for the temperature variance 6 2 is given by (11) with the 
index i replaced by 3. From the algebraic equations for vLu] and u t 6 , Eqs.(9a) and 
(12a), w T e obtain: 
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v9 = —A A 
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Eqs.(15) can be solved using symbolic algebra. The results are: 

, , r (9U dV\ 

{uw ' = "MaT’ a7 


w6 = — Kjj 
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dz 


I<m = ctSm, A'h = er5 w 
5 a/ = — (-s 0 + ^iC// + s 2 Gm) 
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where Gh and Gm are the dimensionless gradients for the mean potential temperature 
and the mean velocity 

Gh = (t A r ) 2 , Gm = ( r 5) 2 
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and 


Z? = do + di G// + d^Gu + dzG 2 H 4- </ 4 G//Gm + d^G 2 ^ 


(18a) 

(186) 

(18c) 
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1 3 

si — — A 4 (A 6 + A 7 ) + 2A 4 A 5 (Ai — -A2 — A 3 ) + -AiA 5 A 8 

3 

S2 — ~Ai(Ag Ay), ^4 = 2A5, S5 = 2A 4 
8 

^6 = “A 5 ( 3 A 3 — Aj) — -Aj A 5 ( 3 A 3 — A2) 4 - — Ai(A6 — A7) ( 19 ) 
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d. Realizability conditions for level 2.5 model 


Realizability conditions are common to all second-order closure models. For the 
present 2.5 level model, there are a few limitations on the values of Sm and Sh , (17a) 
and (17b). Specifically, the two variables Gm and Gh must be limited to certain 
domains outside of which the model may produce unphysical results since some un- 
derlying assumptions (e.g., that departure from isotropy be small) may no longer be 
valid. 

Let us first consider the limitation on buoyancy. Gh may be negative (unstable), 
zero (neutral) or positive (stable). Assuming that production equals dissipation for the 
turbulence kinetic energy e [see Eq. (22) below], and taking the limit Gm —* 0 and 
noticing that Gm is always non-negative, we have 

S h (0,G„)G h + 2 > 0 (20a) 

Substituting Eq.(17b) into Eq.(20a) yields the relation 

— (s 4 + '2d\) HH [(s 4 -f 2oh) 2 — 8dp(s^ + 2df)f^ 

2(^5 + 2df) 

For the model constants used here (see Section 6), this minimum value of Gh is 
— 10.8; the negative value indicates that it occurs in the unstable region. 

Next, we examine the limitation on the shear number. Gm should always be non- 
negative. Following Hassid and Galperin (1983), who argue that an increase of shear 
should not result in a decrease of normalized momentum flux, we apply the following 
condition, 



d ( mv 2 + vw 2 )H 2 

dG m e 


> 0 


Losing Eqs. (16-18), Eq.(21a) can be reduced to a cubic inequality in Gm , 


(21a) 


S2d$G 3 M + [(351^5 — SidfjGn + 3sods — 52 ^ 2 ]^?^ + [(.Sid 4 — ‘5s2df)G 2 H 


+ (s]C?2 sod± — 352^ )Gh — 3^2 do + sodflGM — (^0 ■+■ s\Gn^dzG 2 H ~\~ d\Gn + df) < 0 

( 216 ) 
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Although Eq.(21b) can be solved exactly, one may use the following approximate ex- 
pression based on the fact that the terms containing s? and d 5 are relatively small, 


Gm < 


do + d] Gh + djGjj 
c?2 -f (I4 G 11 


(21c) 


e. Level 2 model 


If we assume that production and dissipation equal each other, the differential 
equation for e, Eq.(14), reduces to 


Sm{Gm , Gh)Gm — Sh{Gm,Gh)Gh — 2 — 0 ( 22 ) 


which can be re-written as an equation for Gm (or for Gh) that depends on only one 
parameter, the gradient Richardson number, 


Ri = 


G_h_ 

G m 


N 2 

!p 


(23 a) 


The resulting equation is: 


(ci-ftf 2 + c^Ri + c 3 )G 2 m + ( c 4 Ri + c 5 )Gm + C6 = 0 (23 b) 


where 

Cj = S $ + 2c? 3, C2 = —5] + 56 + 2d|, C3 = —52 + 2t?5 

C 4 = 54 -f- 2di, C5 = — 5o + 2c?2i C6 — 2do (23c) 

Eq.(23b) is a simple algebraic expression for the variable Gm- Substituting the Gm 
solved from (23b) into (17a-b) we can plot the stability functions Sm and 5// as finction 
of Ri (Figs. 1 and 2). 

The critical Richardson number Ri c , beyond which stable stratification effectively 
suppresses the turbulence, can be found by considering the limit e — ► 0 , i.e., Gm — ► 00. 
In this limit, Eq.(23b) is satisfied only if the coefficient of the quadratic term vanishes, 
which yields 


Ri c 


—£2 + (C2 2 — 4CiC 3 )2 
2cj 


(24a) 


13 


Using the model constants determined in Section 6, we obtain 


Ri c = 0.96 (246) 

Although most previous second order closure models give Ri c ~ 0.2, there is a 
variety of data that are in favor of a Ri c of order one. Early laboratory data by Taylor 
(as cited in Monin and Yaglom, 1971) showed that turbulent exchange exists even when 
Ri > 1. In 1964, Webster’s laboratory measurements showed that mixing persists up 
to Ri ^ 1. In the oceanic PBL, Martin (1985) showed that Ri ~ 1 is needed to obtain 
the correct mixed layer depth at Papa and November stations. More recently, DNS 
(direct numerical simulation, Gerz et ah, 1989) and LES (Wang et al. 1996, Kosovic 
and Curry 2000) show that turbulence exists up to Ri ~ 1. Historically, the criterion 

Ri > - (24c) 

was established by Miles (1961) and Howard (1961) on the basis of linear stability 
analysis. However, when nonlinear interactions were included, Abarbanel et al. (1984) 
showed that the sufficient and necessary condition for stability is not (24c) but that 

Ri > 1 (24c?) 

which is in agreement with our result (24b). 

6. Determination of model constants 

In order to determine the model constants defined in (9b), (10c) and (12c), we 
will employ the related time scale ratio expressions formulated in a recent theoretical 
turbulence model that was based in part on RNG (Renormalization Group) methods 
and whose predictions were tested on different flows (Canuto and Dubovikov, 1996a,b; 
1997): 

y — -El — x i = —A = 0.107, A 5 = — = 5(1 + a to l ) 

r 5 15 r p e 

As = (1 - 7i)— = (1 - 7i)<7to, 7i = 1 (25a) 

T O 
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where cr t0 is the turbulent Prandtl number in neutral flows and will be determined later 
on. To determine A 2 , A 3 and A 4 , we adopt the following expressions (Shih and Shabbir, 
1992; Canuto, 1994): 

qi = 6q 5 , q 2 = -(2 — 7q 5 ), q 5 = — (1 + -F 1 / 2 ) 

F = 0.64, ft = \ (256) 

Substituting (25a) and (25b) in (9b) yields 

A 2 = 0.0032, A 3 = 0.0864, A 4 = 0.1 (25c) 


We parameterize t (the dissipation rate of e) as 

..3 


e = 




(26a) 


which corresponds to 



(266) 


where the dissipation length scale t ~ kz as z ~ 0 and the constant Bi is defined as 
B\ = q 3 /ul where u . is the friction velocity, and the value of B\ must be determined. 
We show that Bi is related to the values of Aj,A 2 and A 3 by assuming a logarithmic 
wind profile near the neutral surface layer (taking the mean wind direction as the x 
direction). The result is: 


=(^A 1 -A 2 + iA 2 )" 3 / 4 = 19.3 


(27) 


To determine the values of A 6 , A 7 and cr t 0 , we need some auxiliary relations. First, 
from (15g,i) an expression for the ratio of the vertical and longitudinal heat fluxes can 
be derived, 


w0 1/2 

~^- XsGm 


T-l 


&t + 2^ 6 


(28a) 


where a t = S M /S H is the turbulent Prandtl number. Webster (1964)’s experimental 
data show that this ratio approaches unity as Ri ~ 0, 


A ,B[ 


- 2/3 


&to + 2^ 6 


-1 


= 1 


(286) 
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where we have used the fact that Gm = B^ 3 at Ri = 0. Second, similar to the 
derivation of (26b), in a near-neutral surface layer, from (15g,i) we obtain 


1 


A 2 — (1 + A2 — 3A 3 )(7toA5 — (^6 ~ A7)(2cr t o + A 6 + A7) — 0 

Using (10c) and (28b-c), we obtain, in the near neutral surface layer, 


1 


4 / 3 / 


(28c) 


4 4 

<23 — 7 + ~&to 
0 0 


l_I 5l 2/3 (l + A 2 -3A 3 ) 


(28 d) 


and A 6 and A7 can be obtained using (28d) in (10c). We still need to determine a value 
for a to in a consistent manner. From the third expression of (25a) and(28b-d), <j t o is 
found to be related to B\, A 2 and A 3 as follows: 


&tQ — 


75 - 3 B 2/3 + 3 1 / 2 

1875 + 150Z?] 2/3 + (403 + 400 \ 2 

- 1200A 3 )5i /3 

1/2 

2B 2/3 

3 + 45j /3 (1 + A 2 -3A ; .)' 



= 0.82 
(28e) 


So it follows that: 


A 5 = 11.04, A 6 = 0.786, A 7 = 0.643, A 8 = 0.547 (28/) 

To summarize, the model constants are determined to be: 

{B x . A,,A 2 .A 3 ,A 4 ) = (19.3,0.107,0.032,0.0864,0.1) 

(A 5 ,A 6 ,A 7 ,A 8 ) = (11.04,0.786,0.643,0.547) (29) 


7. Comparison with Mellor-Yamada model and experimental data 


a. Mellor-Yamada model: a special case 

The MY model (Mellor and Yamada, 1982) corresponds to: 


A! =4 


6(|M 


1 , 3A 1 


A 9 — A 3 — A4 — ~A — 


B x 


\ - B} a — A — 1 \ - T ±- B2 

As-—, A 6 -A 7 -l, A 8 - t _ — 


Thus, 


A = 


6.4 1 

~K' 


o 1 = a 2 = 03 — 71 =0 


(30a) 


(306) 
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where the constants Ai t B\, A 2 and B 2 are determined by Mellor and Yamada to be: 


{A u B u A 2 ,B 2 ) = (0.92,16.6,0.74,10.1) (30c) 

which correspond to a set of value for the model constants in the present model 

(B1,A 1 ,A 2 ,A 3 ,A 4 ) = (16.6,0.168,0.166,0.166,0.166) 

(A 5 ,A 6 ,A 7 ,A 8 ) = (7.48,1,1,0.608) (30 d) 

Substituting (30d) into (19), (23c) and (24a) yields 


Ri c — 0.193 {MY model) 


(30e ) 


One of the deficiencies of the MY model, as Mellor and Yamada pointed out them- 
selves, is that in a neutral surface layer, the model cannot distinguish v 2 and u? 2 , the 
lateral and vertical components of the velocity variance, while experimental data con- 
sistently show that w 2 is always significantly smaller than u 2 . The present model solves 
this problem by incorporating more complete pressure correlations into the model clo- 
sure so that the new model has more freedom to allow v 2 and w 2 to be different. To 
see this more closely, in a neutral surface layer, we reduce (15a-c) of the present model 
to: 


~ — o + 

Q 1 3 


1 A 2 + 3 A 3 


3 

2A 2 

3 


w 2 1 A 2 — 3A 3 
-q 2 = 3 + 3 


(31a) 

(316) 

(31c) 


In the MY model > 2 = A 3 , which makes v 2 = w 2 , while in the present model A 2 and A 3 
are two independent parameters, and we choose to determine them according to Shih 
and Shabbir (1992)’s expressions that are derived from theoretical considerations and 
have been shown to be consistent with measured data. 


b. Comparison with measured and LES data 

The turbulent Prandtl number, a t — I\m/ Kh, is one of the important parameters 
of turbulence. We compare the inverse of o t as a function of the gradient Richardson 
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number Ri resulting from both the present model and the MY model with the experi- 
mental data of Webster (1964). It is clear that turbulence in the stably stratified flow 
exists well beyond the MY critical value Ri % 0.2. According to the experimental data, 
the critical value of Ri should be of order unity, and the present model falls within the 
range of the measured data (Fig. 3). 

We also compare the vertical and lateral heat flux ratio —w6/u6 (as a function of Ri) 
resulting from both the present model and the MY model with the experimental data 
of Webster (1964). Webster described the ratio as "(being) seen to fall catastrophically 
from unity in neutral conditions to only about 0.5 at Ri equal to 0.2 and even less for 
higher Richardson numbers." The present model gives the critical Richardson number 
Ri c = 0.96, in agreement with the data (Fig. 4). 

We then compare the present model with the recent LES of Kosovic and Curry 
(2000) on a stably stratified PBL, at hour 12 of the high-resolution case NLHRB, 
when a quasi-steady state is reached. In our simulation we use the level 2 model 
since we are particularly interested in the behavior of the model when the gradient 
Richardson number Ri varies; for the length scale formula we use 


*o = 0.1 


Jo” zqdz 
Jo” qdz 


Ci = 


kzCq 

Cq + KZ 


t=l 


Ci 


§f < 0 


min{l\ , 0. 53£) : g>0 


(32) 


One can see from figures 5-7 that the present model is more consistent with the LES 
result than the MY model, as we will discuss below. 

It is very informative to examine the nondimensional shear and potential temper- 
ature gradients defined as 


k z _ 

<f> — — s 

^ m — *'■'5 

u w 


KZU , <90 

\w0s\ dz 


(33) 


where u, and w6 s are the friction velocity and the surface potential temperature flux 
respectively and 5 is the shear given by Eq.(18b). In Figs. 5 and 6 we plot l/$ m and 
l/$h as function of Ri. The graphs indicate that the present model can reproduce 
the observed Kansas data as analyzed by Businger et al. 1971) in the context of 
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Monin-Obukhov similarity theory for Ri < 0.2. For Ri > 0.2, the present model shows 
that turbulence still exists although is weaker, in agreement with Kosovic and Curry 
(2000)’s LES data. This result is also in agreement with Brown et al. (1994) and 
Andren (1995). In the same figures we also plot the results of the MY model, which 
deviate from the Kansas data for Ri < 0.2, and fail to reproduce the turbulence beyond 
Ri — 0.2 found in the LES. 

The PBL height is one of the most important quantities in any PBL modeling. The 
PBL height is usually defined as the height at which the turbulent kinetic energy or the 
magnitude of the momentum flux decreases to a small fraction of the corresponding 
surface value ; or it may be defined as the height at which the (positive) temperature 
gradient reacheas a certain value from below. In any case, the top of the PBL lies in 
a region where the turbulence is stably stratified and, given the mean profiles of the 
wind and the temperature (and thus given Ri), a higher intensity level of turbulence 
yields a greater PBL height. The MY model, however, underestimates the PBL height 
(Yam ad a and Mellor 1975). Since the present model predicts larger critical Richardson 
number and produces more turbulence for a given Richardson number, greater PBL 
heights can be achieved (Fig. 7). 

8. Conclusions 

With a single new input, the most updated expressions for the pressure- velocity and 
pressure-temperature correlations, we have derived a second-order closure turbulence 
model to describe the PBL. One of the main features of the new model is that it yields a 
critical Richardson number of order unity, rather than ~ 0.2 as given by most previous 
models. The larger critical Richardson number is in agreement with measured and LES 
data and the stability analysis that includes non-linear interactions. The new model 
reproduces the Kansas data as analyzed by Businger et al. (1971) for Richardson 
numbers smaller than 0.2 as well as the LES and laboratory data for Richardson 
numbers up to unity. Another improvement is that the present model allows a match of 
the surface data, which was not possible in previous second order closure PBL models. 
In addition, the new model produces greater PBL height than the previous models. 
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Figure Captions 


Fig.l. The stability function Sm versus the gradient Richardson number Ri. The 
solid line represents the present model; the dotted line, the MY model. 

Fig. 2. The stability function 5// versus the gradient Richardson number Ri. The 
solid line represents the present model; the dotted line, the MY model. 

Fig. 3. The inverse turbulent Prandtl number <j t 1 (normalized by its value for 
neutral stratification) versus the gradient Richardson number. The solid line is the 
result of the present model at level 2. The dotted line represents the level 2 MY 
model. The experimental data by Webster (1964) are redrawn here as filled circles. 
The present model yields a much larger critical Richardson number (ss 1) than the 
Mellor-Yamada model (ss 0.2). 

Fig. 4. Ratio of the rates of heat transport in the w-direction (vertical) and the 
u-direction (horizontal, along the mean flow), —w9/u6, versus the Richardson number. 
The solid line represents the result of the present model, while the dotted line represents 
the MY model. The experimental data (Webster, 1964) are redrawn here as filled 
circles. 

Fig. 5. N ondimensional shear as a function of the gradient Richardson number. The 
crosses represent the LES simulation of Kosovic and Curry (2000), case NLHRB (their 
Eq. 27). The solid line represents simulation results using the present model, while 
the dotted line, simulation results using the MY model. The dashed line represents 
the measured Kansas data (Businger et ah 1971). 

Fig. 6. Similar to Fig. 5 but for the nondimensional potential temperature gradient 
(Eq. 28 of Kosovic and Curry, 2000). 

Fig. 7. PBL height as a function of the dimensionless time tf c , where f c is the 
Coriolis parameter. Cross: LES result; solid line: present model result; dotted line: 
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MY model result. 


Appendix 

The Level 3 PBL model 


In the level 3 PBL model, the turbulent temperature variance 6 2 is solved from its 
prognostic equation (instead of from an algebraic equation): 


d 


oe. 


e 2 


— 0 2 + —w 6 2 = -2 —w 6 - 2— 
Dt dz dz to 

From (10a), the algebraic equation for the heat flux w0 is: 


(AF) 


w 6 — — Ar 1 ' 


+ i(Ae - A 7 ) + ^voj 


+ A 5 1 A 0 gctrO 2 


(A 2) 


All the other algebraic equations for the Reynolds stress and the heat flux are the same 
as the level 2.5 model, (15a*h), except (15i), which is replaced by (A2) in the level 3 
model. We solve (15a-h) and (A2) using symbolic algebra and the results are: 


(uuqme) = -erS M I — , ^ 


dU dV\ 


„ „ oe 

W0 = —€T Stf—- + 7 C 
dz 


(A3) 


where 


lc = 


3[A 5 + A 4 G„ + A 5 (Ag - LA l)G M ] 

D 


\ o gQT0 2 


(.44) 


is the counter-gradient term which is absent in the level 2.5 model and D is the same 
as in (18c). The structure of the stability function Sm differs from the Sm in the level 
2.5 model (17a) by an extra term, 




■So + -Sj G h + + S3\o(gaT) 


, 9 2 


(A 5) 


where 

S 3 = -A 4 [2A 5 (A2 + 3A3) -f 3(A6 + A7)] (-46) 

The form of the stability function Sh remains the same as that of the level 2.5 model 
(17b). The expressions for the model constants are the same as in (19) except that 
in level 3 model A s = 0. In addition, the constant A 0 needed in (A4) is set to 2/3 
according to (10c) and (25a). 
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